What does measure the scaling exponent of the correlation sum in the case of human 

heart rate? 
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It is shown that in the case of human heart rate, the scaling behaviour of the correlation sum 
(calculated by the Grassberger-Procaccia algorithm) is a result of the interplay of various factors: 
finite resolution of the apparatus (finite-size effects), a wide dynamic range of mean heart rate, the 
amplitude of short-time variability being a decreasing function of the mean heart rate. The value of 
the scaling exponent depends on all these factors and is a certain measure of short-time variability 
of the signal. 

PACS numbers: PACS numbers: 87.10.+e, 87.80.Tq, 05.45.-a, 05.45.Tp 



Heart rate variability (HRV) is often thought to be 
driven by deterministic chaos inside the heart. This chaos 
is explained in terms of dynamics of the complex of sino- 
atrial and atrio-ventricular nodes, which has been suc- 
cessfully modeled as a system of non-linear coupled os- 
cillators, responsible for the heart rhythm 0,0. 

As a consequence, correlation dimension and related 
quantities (like Lyapunov exponents and Kolmogorov en- 
tropy, etc.) have been thought to be important non-linear 
measures of HRV. In particular, Babloyantz and Des- 
texhe concluded Q that high values of the correlation 
dimension indicate the healthiness of the heart. Later, 
the correlation dimension of the heart rate signal has 
been calculated in a vast number of papers. Meanwhile, 
it has been pointed out that physiological time-series are 
typically non-stationary and noisy, and therefore, the cor- 
relation dimension cannot be calculated reliably 0, IE 0| 
- the fact which is now widely accepted. In the case 
of human heart, the "noise" comes from the autonomous 
nervous system in the form of inputs regulating the heart 
rate (cf. 0,1111). These mostly non-deterministic signals 
suppress almost completely the underlying determinis- 
tic signal. Futhermore, it has been emphasized that a 
reasonable fitting of a correlation sum to a power law 
does not necessarily mean that the obtained exponent 
is the correlation dimension of the underlying dynami- 
cal system; instead, thorough non-automatable verifica- 
tion procedure has to be done [ljj. All this leads us to 
the conclusion that the formally calculated correlation 
dimension of a heart rhythm does not correspond to the 
dimensionality of an intrinsic attractor. 

Thus, there are two important observations: (a) the 
correlation sums of human heart rate follow typically a 
scaling law; (b) in most cases, the scaling exponents are 
not the correlation dimensions. Then, a natural question 
arises, what is the physical meaning of these formally 
calculated exponents? 



Our answer to the posed question is based on very 
simple observations, which are valid for healthy patients: 
(a) the long-time variability of the inter-beat intervals 
(around 500 ms) is typically much higher than the vari- 
ability on the time scale of few heart beats (less than 50 
ms); (b) for those periods, when the mean heart rate is 
high (i.e. when the subject is performing physical exer- 
cise) the heart rate variability is low; (c) the heart rate 
is controlled by non-deterministic and effectively random 
signals arriving from the autonomous nervous system. As 
a consequence, in time delay coordinates, an HRV time- 
series generates a baseball bat-shaped cloud of points. 
Although the theoretical value of the correlation dimen- 
sion of such a cloud is infinite, the finite resolution of 
the recording apparatus, finite length of the time-series, 
and the linear structure of the cloud result in a smaller 
value. This is evident for a very narrow "bat" , which is 
efficiently one-dimensional. In what follows we show that 
the correlation dimension reflects the geometrical size of 
the cloud of points. 

The layout of the paper is as follows. First, we give the 
details of the HRV database used for this study. Second, 
we provide a short overview of the research results related 
to the correlation dimension of human heart rhythm. 
Third, we construct a simple model-time-series, the cor- 
relation sum of which scales almost identically to that of 
real HRV data. Finally, we discuss the universality and 
implications of our model. 

The experimental data analyzed in this paper have 
been recorded at Tallinn Diagnostic Center. The record- 
ings of ambulatory Holter-monitoring (24 hours, approxi- 
mately 100 000 data points) were obtained during regular 
diagnostical examinations and covered over 200 patients 
with various clinically documented diagnoses (including 
also many healthy patients). The main groups of patients 
are shown in Table 1 . The resolving power of recordings 
was 6 ms (sampling rate of 180 Hz). The diagnostics and 
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data verification has been made by qualified cardiologist; 
the data preprocessing included also filtering of artifacts 
and arrhythmias. 





Healthy 


IHD 


SND 


VES 


PCI 


RR 


FSK 


No. of patients 


103 


8 


11 


16 


7 


11 


6 


Mean age 


45.5 


65.4 


50.0 


55.9 


47.3 


55.5 


11.7 


Std. dev. of age 


20.5 


11.4 


19.3 


14.3 


11.6 


14.4 


4.6 



TABLE I: Test groups of patients. Abbreviations are as 
follows: IHD - Ischemic Heart Disease (Stenocardia); SND - 
Sinus Node Disease; VES - Ventricular Extrasystole; PCI - 
Post Cardiac Infarction; RR - Blood Pressure Disease; FSK - 
Functional Disease of Sinus Node. 

The concept of correlation dimension, introduced by 
Grassberger and Procaccia [Tj , is designed to reflect the 
number of degrees of freedom of a deterministic system 
(or the dimensionality of an attractor, which, in principle, 
can be fractal). For empirical time-series, the phase vari- 
ables are typically not known. It is expected that the at- 
tractors in the phase space are topologically equivalent to 
the attractors in a reconstructed phase space with time- 
lag coordinates {x(t) 7 x(t + r), . . . x[t + (m — 1)t]}, as long 
as the embedding dimensionality m (the dimensionality 
of the reconstructed phase space) exceeds the dimension- 
ality of the attractor D; here x(t) is the signal, and r is 
a reasonably chosen time lag. This circumstance is ex- 
ploited by the Grassberger-Procaccia method ^lj f° r the 
calculation of the correlation dimension. To begin with, 
the second order correlation sum is defined as 



C 2 (r) = 



N(N - 1) 



5>(r- 



I), 



(1) 



where 9(r) is the Heaviside function, and r, = 
{x(ti), x(ti + r), . . . , x[t + (m — 1)t]}, is a point in the 
reconstructed phase space, and i = 1,2, .... N counts the 
moments of discretized time. For small r, the correlation 
sum is expected to scale as Ci (r) oc r° 2 , assuming that 
Z?2 < m- The exponent D = D2 is called the correlation 
dimension of the system. 

A non-linear dynamical system may be chaotic and 
then the phase trajectory fills the entire phase space. 
In that case, the correlation dimension D2 is equal to 
the number of degrees of freedom (the dimensionality of 
the phase space minus the number of conservation laws). 
This is why D2 is often considered as a measure of the 
complexity of the system. Babloyantz and Destexhe Q 
studied the correlation dimension of the sequence of NN- 
intervals (intervals between normal heartbeats) of human 
heart rhythm. For healthy patients and data series con- 
sisting of 1000 intervals, they found D = 5.9 ± 0.4. It is 
widely recognized that life threatening heart pathologies 
lead to the reduction of the complexity of the HRV sig- 
nal, c.f. [l^. Correspondingly, the correlation dimension 
of the heart rate has been often considered as a measure 
for the healthiness of the heart. 

However, the heart is not an isolated system. Although 
the heart rhythm is generated by the complex of oscilla- 



tory elements, its rate is controlled by non- deterministic 
inputs arriving from the autonomous nervous system. In 
particular, these inputs lead to the increase of heart rate 
when the subject is under physical stress, and to slowing 
down when the subject is at rest, c.f. Healthy heart 
responds easily to these signals, and is able to adapt to 
a wide range of beating rates. This responsiveness gives 
rise to the high variability of the heart rate. Severe heart 
diseases decrease the responsiveness of the heart with re- 
spect to the whole spectrum of signals arriving from the 
autonomous nervous system; this leads to the loss of ap- 
parent complexity of the HRV signal. 

The heart is more responsive with respect to the signals 
of the autonomous nervous system when the heart rate 
is slow, i.e. when the patient is at rest. In that case, the 
heart rate variability is driven by weaker signals, like the 
ones generated by respiration and blood-pressure oscilla- 
tions. These two stimuli are quasi-periodic, the periods 
being respectively a few and 10-20 seconds. It should be 
noted that respiration can be mode-locked to the heart 
rate. This mode-locking has been demonstrated by si- 
multaneous recording of ECG and respiration activity, 
together with the technique called cardiorespiratory syn- 
chrogram 0] . The ratio of the mode- locked periods can 
be small, 2:1, 3:1, 5:2, etc., and the phenomenon can give 
rise to certain patterns in the reconstructed phase space. 
These patterns can be easily misinterpreted as traces of 
an attractor of a non-linear deterministic system, there- 
fore we discuss this aspect in more details. 

As mentioned above, HRV signals generate a baseball- 
bat-shaped cloud of points in a time-lag space. For cer- 
tain patients, the presence of less densely populated satel- 
lite clouds can be observed, see Fig. 1 (note that due to 
the sparse population of these clouds, the presence or 
absence of the clouds has almost no effect on the value 
of D). We were able to show that for all such patients 
of our database, the satellite clouds were not related to 
a deterministic chaos inside the heart. To this end, we 
studied the fluctuation function 



F(u) = (\t„-t n+v 



(2) 



(angular braces denote averaging over n). Unlike in the 
case of "single-cloud-patients", the fluctuation function 
of the patients with satellite clouds revealed a presence 
of an oscillatory component, see Fig. 2a. By dividing the 
entire 24-hour HRV record into one-hour intervals, and 
calculating the amplitude of the oscillatory component of 
the fluctuation function for each interval, we were able 
to locate the periods responsible for the satellite clouds 
in the reconstructed phase space, see Fig. 2b. These 
were always the periods before falling asleep, around 10 
or 11 pm, characterized by a low heart rate and a high 
respiration-driven short-time variability. The phase be- 
tween the heart rate and respiration is locked during tens 
of seconds, confirming the observations of Kurths et al. 
|13|| . Thus, in a certain sense, the heart and respira tory 
complex act as a system of coupled oscillators c.f. [l4j : 
however, there is no determistic chaos inside the heart. 
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Note that our method of mode-locking detection is very 
simple, does not require synchronous respiration rhythm 
recording (unlike the thourough method 0), and can 
be conviniently used to find relatively short (> lOmin) 
locking periods from a 24-hour recording. 
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FIG. 1: Two-dimensional intersection of 3-dimensional recon- 
structed phase space for a patient with pronounced mode- 
locking between heart rate and respiration. The number of 
points per unit cell is given in gray-scale coding. 
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FIG. 2: Patient with 3:1 mode locking between heart rate and 
respiration: (a) heart beat intervals (in milliseconds) plotted 
versus the beat number. Heart rate has a pronounced oscilla- 
tory component; vertical lines mark the period of three heart 
beats, horizontal lines indicate the sequences with coherent 
phase, (b) Fluctuation function (arbitrary units) is plotted 
versus the time lag v (in heart beats); the oscillating compo- 
nent is magnified. 

Our model of the heart rhythm generation is as fol- 
lows. The non-linear deterministic dynamics inside the 
heart is almost completely suppressed by the signals ar- 
riving from the autonomous nervous system. These sig- 
nals control the mean heart rate, but obey also a noise- 
like component, the amplitude of which decreases with 
increasing mean heart rate. In the case of correlation 
dimension analysis, this noise- like component is indistin- 
guishable from a Gaussian noise. Therefore, theoreti- 



cally, the correlation dimension is infinite. The reported 
relatively small values of the correlation dimension are to 
be attributed to the finite length of the time series and, 
most importantly, by finite resolution of the recording 
apparatus. Regarding the length of the record: typically, 
the correlation dimension has been found to be at the 
limit (or beyond) of a credible analysis 0, ■ Indeed, 
it has been suggested jj, [15J that the calculation of the 
correlation dimension D is reliable, if the number N of 
data-points in the time series 



N > 10 D / 2+1 



(3) 



In Table 2, this criterion is compared with the data of 
some papers. 





Ref. [3] 


Ref. [5] 


Ref. [16] 


Ref. [17] 


Correlation dimension 


5.5 6.3 


9.6 10.2 


2.8-5.8 


4-7 


Length of the data set 


10 a 


2- 10 4 


10 4 


2 ■ 10 4 


Required length 


10 4 


10 b 


10 4 


3- 10 4 



TABLE II: Data from papers devoted to the correlation di- 
mension analysis: experimental values of correlation dimen- 
sion, lengths of the underlying data sets, and data-set lengths 
required by Eq. 3. 

In order to test our hypothesis we aimed to construct 
such random time series (using an algorithm as simple as 
possible), the correlation sum of which is similar to the 
correlation sums of the time series of real patients. First 
we analyzed the sequences of NN- intervals extracted from 
ECG recordings. The scaling exponent was calculated 
according to the Grassberger-Procaccia algorithm. The 
six-dimensional embedding phase space was used for cal- 
culations. The choice of the embedding dimensionality 
was motivated as follows. To begin with, the analysis 
for phase space with m > 6 is not reliable due to sparse 
clouds of points (see Eq. 3). Further, m — 6 does still 
make sense, because most of the previous studies have 
reported D < 6. 

Reliable correlation sum analysis is possible only for 
more or less stationary time series, cf. Meanwhile, 
HRV signal is highly non-stationary, as is evidenced by 
the multifractal structure of its long-time dynamics [18j . 
The most stationary period in the heart rate dynamics 
is the sleeping time. This is why we studied only the 
nocturnal part of the HRV records. The scaling exponent 
was determined as the slope of the correlation sum C2 (r) 
in log-log plot by performing root-mean-square fit for the 
almost linear part (at small values of r) of the curve, see 
Fig. 3. Note that the leftmost horizontal part of the 
curve is due to the limited resolving power (6 ms) of the 
medical equipment: if two NN- intervals differ less than 
6 ms, they are recorded to be of the same length. The 
scaling exponents ranged from D — 4.2 to D = 5.1 and 
were almost uncorrelated with the diagnoses (see Table 

Further we generated two surrogate time-series with 
Gaussian noise: (i) plain Gaussian noise added to a mean 
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p,% 


IHD 


SND 


VES 


PCI 


RR 


FSK 


Healthy 


89.4 


21.9 


3.5 


18.4 


2.4 


71.5 


IHD 




34.1 


12.0 


17.6 


7.1 


69.4 


SND 






66.8 


52.9 


45.7 


54.4 


VES 








73.0 


67.6 


25.2 


PCI 










95.7 


26.7 


RR 












15.9 



TABLE III: p-values of the Student test for the seven group 
of patients. Abbreviations are explained in Table 1. 

rate (see Fig. 4c); (ii) time-series with variable mean rate 
and modulated noise, generated according to the formula 

t n = a + 6sin(/n) + cg(n)y/l.l + sin(/ra), (4) 

see Fig. 4b. Here, t n denotes the duration of n-th in- 
terval; g(n) is a random normally distributed value with 
zero mean and standard deviation equal to 1. The term 
6sin(/n) models the variability of the mean heart rate 
due to physiological processes (physical activity, blood 
pressure oscillations, etc.). The term ylT + sm (f n ) re- 
flects the empirical observation that the short-time vari- 
ability of heart rhythm increases together with the mean 
heart beat interval. Note that both the square-root and 
sinusoidal dependances are rather arbitrary, the model is 
not sensitive neither with respect to the particular func- 
tional dependencies nor with respect to the modulation 
frequency /. The numerical values of these parameters 
have been chosen as follows: a = 500 ms, b = 110 ms, 
/ = 0.005, c = 3.5 ms; the values of t n were rounded to 
the nearest multiple of 6 ms (the "resolving power" ) . 

For a Gaussian signal, the correlation dimension is in- 
finite and the scaling exponent should be equal to the 
embedding dimension m = 6. This is exactly what is ob- 
served for plain unmodulated Gaussian time-series, see 
Fig. 5. However, for the noise of modulated amplitude, 
the finite size effects are significant. Depending on the 
parameters b and c, the correlation sum in Fig. 6 can be 
almost indistinguishable from the ones of real patients, 
see Fig. 4. The scaling exponent D of such time-series 




FIG. 3: Heart beat intervals (in arbitrary units) are plotted 
versus the beat number: (a) a real patient; (b) surrogate data 
(modulated Gaussian noise); (c) plain Gaussian noise added 
to a constant "heart" rate. 




FIG. 4: Correlation sums of a healthy patient, a plain Gaus- 
sian signal, and a modulated Gaussian signal in logarithmic 
plot. Embedding dimensionality m — 6. 



depends on the amplitude of the Gaussian noise and on 
the resolving power (which, unlike in the case of real ap- 
paratus, was also freely adjustable). By adjusting the 
above defined parameters b, c, and the resolving power, 
we were able to obtain the values ranging from D = 4 to 
L> = 6. 

In conclusion, comparative analysis of real and surro- 
gate data confirmed that the scaling behaviour of the cor- 
relation sum and finite values of the formally calculated 
correlation dimension D are the result of the interplay 
of the following factors: finite resolution of the recording 
equipment (which leads to finite-size effects), a signifi- 
cant level of long-time variability (the dynamical range 
of the mean heart rate exceeds the typical level of short- 
time variability), the fact that the amplitude of short- 
time variability is a decreasing function of the mean heart 
rate. The value of the scaling exponent is mostly defined 
by the dynamics of the short-time variability and resolv- 
ing power of the recording apparatus; it depends also on 
on the maximal embedding dimensionality. Therefore, it 
can be used as a certain measure of short-time variabil- 
ity of the signal (however, in order to obtain comparable 
values, time-resolution and record length have to be kept 
constant). The diagnostic and/or prognostic value of this 
measure is possible, but found to be non-significant for 
our patient groups (see Table 2). We have also shown 
that the above drawn conclusion remains valid even in 
these cases, when certain patterns (satellite clouds) can 
be observed in time delay coordinates of heart rhythm 
(see Fig. 3). These patterns are typically due to the res- 
piratory sinus arrhythmia and mode coupling between 
respiration and heart rhythm without any relation to the 
(possibly) derministic dynamics inside the heart. Finally, 
we have devised a simple method of detecting the pres- 
ence of this mode coupling, based on fluctuation function 
(2). 
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